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Abstract. In this paper, we present a stabilized mixed formulation for unsteady Brinkman equa- 
tion. The formulation is systematically derived based on the variational multiscale formalism and 
the method of horizontal lines. The derivation does not need the assumption that the fine-scale 
variables do not depend on the time, which is the case with the conventional derivation of mul- 
tiscale stabilized formulations for transient mixed problems. An expression for the stabilization 
parameter is obtained in terms of a bubble function, and appropriate bubble functions for various 
finite elements are also presented. Under the proposed formulation, equal-order interpolation for 
C^) ' the velocity and pressure (which is computationally the most convenient) is stable. Representative 

numerical results are presented to illustrate the performance of the proposed formulation. Spatial 
' and temporal convergence studies are also performed, and the proposed formulation performed well. 

1. INTRODUCTION 

(N 

> ; Among the engineering and environmental issues that demand attention at present, enhanced oil 

recovery (EOR), carbon-dioxide sequestration, and contaminant transport in heterogeneous media 

00 ■ 

£f) • can lay claim to pride of place. The models currently used are mostly steady state models and 

reliable computational techniques do exist for them. However, there is a big drawback because it 
is only if the transient effects are considered that it will be possible to later study the coupling 
phenomena that occur due to interaction between different media. Our intention is to propose a 
stabilized mixed formulation for an unsteady model, since the existing techniques have their own 
disadvantages. However, before developing the finite element formulation, one needs to turn to 
theory of mixtures to trace the development of the models. 

The foundations of the theory of interacting continua (also known as theory of mixtures) were 
laid down by Truesdell [U [2] . One particular physical situation in which it can be gainfully employed 
is encountered frequently in the flow of a fluid through a porous solid. Ideally, one would like to 
be able to predict all quantities of interest in such a situation knowing the nature of the mixture 
constituents (solid and fluid in this case). But this is easier said than done, for the deformation 
of the solid due to the flow of fluid through it could be highly non-linear, and the fluid itself may 
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show non-linear behaviour. Being able to accurately describe such situations quantitatively is the 
ultimate goal of theory of mixtures and it has been shown in [3] that starting from the balance 
of linear momentum, one can derive a whole gamut of mathematical models pertinent to different 
situations with varying levels of generality. However, using these models to actually solve realistic 
engineering problems is another matter altogether. Developing computational techniques to address 
such models is one of the important thrusts of current research. Herein, we take into account the 
transient and viscous effects of the fluid (while still neglecting the motion of the solid) , which gives 
rise to unsteady Brinkman equation. The main intent of the present article is to present a stabilized 
mixed formulation for the solution of the unsteady Brinkman equation. 

Darcy's equation [3J describes the flow of a fluid through a porous solid due to pressure gradients 
when a wide variety of assumptions are justified. It has been used to study various phenomena 
arising in groundwater hydrology [5j[6], enhanced oil recovery 0|8], carbon-dioxide sequestration 
[§|I1U] to name a few. However, it should be noted that Darcy's equation is merely an approximation 
to the balance of linear momentum for a fluid as it flows through a rigid porous solid (see References 
[31 [TTj [T21 [T5] for a detailed derivation using the theory of interacting continua, which is also known 
as theory of mixtures). 

Rajagopal has shown that, within the context of theory of mixtures, the Darcy equation can 
be obtained after making a plethora of assumptions about the solid and the fluid. The assumptions 
include neglecting the deformations of the solid, assuming the flow to be steady, assuming a special 
form for the drag due to the friction at the pores, and many more. The restrictions under which the 
Darcy equation was derived are rather stringent, and in the assumptions are systematically 
relaxed and an hierarchy of models is thus obtained, where the Darcy model is the one derived under 
the most restrictive assumptions. Moreover, it has been demonstrated in References [141 IS] that 
for a class of flows involving high pressure gradients, the Darcy equation is a poor approximation. 
However, despite its drawbacks, it remains the most popular model to describe the flow of fluids 
through porous solids, and for a large class of flows, it remains applicable. In the remainder of the 
paper, unless otherwise mentioned, the Darcy equation tacitly refers to the steady state equation. 

If the viscous effects in the fluid are deemed important, as they are in many cases, and taken 
into account and modeled as those of the familiar Navier-Stokes fluid, the eponymous Brinkman 
model is obtained [16]. By further relaxing the assumption in the Brinkman model that the flow be 
steady, we take transient effects into account while still neglecting the inertial non-linearities, and 
one thus obtains (as in [11]) the unsteady Brinkman equation, the solution of which is the subject 
of the present article. 

If one desires to use the finite element method to solve the Brinkman equation (or Darcy equa- 
tion), then the celebrated Ladyzhenskaya-Babuska-Brezzi (LBB) condition |17t I18| must be either 

2 



satisfied or circumvented. The LBB condition imposes severe restrictions on the classical mixed 
formulation with respect to the order of interpolation for the pressure and velocity in problems 
involving incompressibility. In particular, the classical mixed formulation is unstable for the equal- 
order interpolation for the velocity and pressure, which is computationally the most convenient. 
There is thus a necessity for stabilized mixed formulations, especially those that are stable under 
the equal-order interpolation for the velocity and pressure. Variational Multiscale (VMS) formalism 
provides a systematic way of developing stabilized formulation |19j . Some of the earlier works on 
mixed methods applied to flows through porous media are [20 \ \21 \ [22l 1231 124 j 25, 26J. A thorough 
discussion of stabilized methods is beyond the scope of this paper, and a reader interested in this 
subject should refer to [271 EH [SH [30l EQl EH [321 [33] . 

Remark 1.1. While it is possible to avoid a mixed formulation by using primal or single field 
formulations in case of the Darcy equation, such a technique cannot be applied to the unsteady 
Brinkman equation. Moreover, the primal formulation has the disadvantage of poor approximation 
of velocity for low-order finite elements (see [3~4"1 13"5] ). 

This VMS method for both Brinkman (or Darcy's equation) involves decomposing the velocity 
field into coarse/resolved scales and fine/unresolved scales. Modeling of the unresolved scales leads 
to a multiscale/stabilized form of the corresponding equation (see Reference [36J for details). 

The subsequent attempts, quite naturally, were to use the VMS method to solve the correspond- 
ing unsteady problems. The traditional way of treating the time derivative in the unsteady equation 
(unsteady Darcy/Stokes/Navier-Stokes) has been to assume that the fine/unresolved scales are in- 
dependent of time while allowing the coarse/resolved scales to depend on both the spatial and 
temporal variables |37[ I38|. However, such an assumption is philosophically undesirable as it seems 
motivated by the need to somehow solve the fine scale problem rather than a sound physical basis. 
For completeness, we have outlined this traditional way (using a semi-discrete method) of deriving 
stabilized formulation in Appendix 17.21 

Another technique that has been used to treat such unsteady problems is the so called "space- 
time finite elements," where the temporal domain is also discretized via shape functions just like 
the spatial domain |39|, [4"0l I41j . However, the disadvantage with this method is that one needs a 
four-dimensional mesh for a (spatially) three-dimensional problem, the implementation is not as 
convenient, and there are also the same issues with post-processing when visualizing the numerical 
results of a (spatially) three-dimensional problem. 

The method followed in this article sidesteps both the preceding issues discussed in [371 EE] and 
[321 SOI H] • We use Rothe's method [42J (also called the method of horizontal lines) and avoid the 
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mathematical legerdemain involved in justifying the assumption of the time independence of the 
fine/unresolved scales. 

The terminology "method of lines" when used with reference to numerical techniques is usually 
understood to imply method of "vertical" lines as distinct from the method of "horizontal" lines 
used in this paper. The method of lines |43t |4"4"] is a favoured numerical technique for the solution 
of parabolic partial differential equations (PDEs) wherein the PDE is discretized in all but one 
dimension and thus reduced to a system of ordinary differential equations (ODEs). On the other 
hand, in the method of horizontal lines, the temporal derivative is discretized first (usually by a 
finite difference scheme) and the original parabolic PDE is thus reduced to a (time-discretized) 
elliptic PDE. The entire machinery of techniques and tools used to study elliptic PDEs (which are 
quite well understood) thus becomes available and therein lies the advantage of the method [45J. 
In a different context, the method of horizontal lines has been used before. In |46j for instance, a 
finite element formulation for transient incompressible viscous flows stabilized by local time-steps is 
derived. However, the central idea of the present paper is different; to use the method of horizontal 
lines in conjunction with the VMS formalism and derive a consistent formulation which circumvents 
an ad hoc assumption that is currently made. 

1.1. Main contributions of this paper. In this paper, we systematically derive a stabilized 
mixed formulation for unsteady Brinkman equation using the method of horizontal lines and the 
variational multiscale formalism. The stabilization terms and parameter are obtained in a math- 
ematically consistent manner. An important feature is that the derivation does not need the 
assumption that the fine-scale variables do not depend on the time, which is the case with the 
conventional derivation of multiscale stabilized formulations for transient mixed problems. Under 
the proposed formulation, equal-order interpolation for the velocity and pressure (which is com- 
putationally the most convenient) is stable. We also derive a stabilized mixed formulation for the 
unsteady Darcy equation as a special case. We also show through numerical experiments that the 
proposed stabilized formulation posses good spatial and temporal convergence properties. 

1.2. An outline of the paper. In Section[2l we outline the governing equations. In Section[3l we 
present a stabilized mixed formulation based on the variational multiscale formalism and the method 
of horizontal lines for the unsteady Brinkman equation. In Section UJ the stabilized formulation 
obtained earlier is specialized for the case of the unsteady Darcy equation. Representative numerical 
results are presented in Section and conclusions are drawn in Section 
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2. GOVERNING EQUATIONS: UNSTEADY BRINKMAN MODEL 

Let f2 C M. nd be an open bounded domain, where u nd" is the number of spatial dimensions. Let 
r denote the boundary, which is assumed to be smooth. A position vector is denoted by x G Cl, 
where an over-bar denotes the set closure. The time is denoted by t € [0, T], where T is the total 
time of interest. The gradient and divergence operators with respect to x are, respectively, denoted 
by grad[-] and div[-]. The velocity (vector) field is denoted as v(x,t), and the pressure (scalar) field 
is denoted as p(x,t). The coefficient of viscosity, drag, and density are respectively denoted by p, 
a, and p. The prescribed specific body force is denoted by b(x,t), while the prescribed velocity on 
the boundary is denoted by v p (x, t). The unsteady governing equations can be written as follows: 

dv 

p— = — gradfc] — av + divht gradk?]] + pb in Q x (0,T) (la) 
at 

div[u] = infix(0,T) (lb) 

v(x,t) = v p (x,t) onIx(0,T) (lc) 

v(x,0) = vq(x) in Q (Id) 



where vq(x) is the prescribed initial velocity. Few remarks are in order about the mathematical 
model given by equations (fTa|) - (fTd|) . 

Remark 2.1. The drag coefficient is related to the coefficient of viscosity through a = — , where 

rZ 

k is the coefficient of permeability. One obtains unsteady Stokes ' equation if the drag term av is 
neglected. On the other hand, one obtains unsteady Darcy equation if the term div[p grad[i>]] is 
neglected and only the normal component of the velocity on the boundary is prescribed. 

It is noteworthy that a more general form of transient Darcy- Brinkman equation can be written 
as follows: 

p— = — grad[p] — av + div[^ grad[t>]] + pb (2) 

where the factor p can be different from the coefficient of viscosity. Physically, the friction at the 
pores of the solid as the fluid flows past it is modeled as a(v — i> so iid)j which reduces to av under 
the assumption that the solid is rigid and does not undergo any motion. On the other hand, the 
viscous effects within the fluid itself are modeled as div[p grad[i>]]. There are adequate discussions 
in the literature (including the paper by Brinkman |47j ) about how the factor p is related to the 
coefficient of viscosity p. However, for a medium with high porosity, it is reasonable to assume that 

p=pM- 
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Remark 2.2. It is important to note that another constraint must be enforced in order to have a 
unique solution for the pressure. This is usually enforced by demanding that 

P dQ = (3) 

Jn 

In numerical simulations, however, it is more convenient to prescribe the pressure at a point, which 
is employed in this paper. 

Remark 2.3. The balance of linear momentum for the fluid, on neglecting the non-linear convective 
term, takes the form 

p^ = div[T]+pb (4) 

Appealing to frame indifference (or even Galilean invariance), the partial stress of the fluid T is 
actually modeled as 

T = -pi + fj, (gra<%] + grad[u] T ) (5) 

It is only when fx is a constant that the incompressibility constraint div[v] = allows the balance 
of linear momentum to be simplified further and it takes the form as given in equation (jla[) . This 
point becomes crucial in developing (realistic) nonlinear models in which the viscosity depends on 
the pressure. 

The main aim of this paper is to develop a stabilized mixed formulation for the above out- 
lined system of partial differential equations (|la|) — (|ldj) . It is well-known that the semi-discrete 
formulation based on the classical mixed formulation has numerical deficiencies. In particular, the 
equal-order interpolation for the velocity and pressure (which is computationally the most conve- 
nient) is not stable. In the next section we develop a stabilized mixed formulation for unsteady 
Brinkman equation under which the equal-order interpolation for velocity and pressure is stable. 
An important aspect is that the derivation does not need the assumption that fine-scale variables do 
not depend on the time, which is the case with the conventional derivation of multiscale stabilized 
formulations. 

3. A STABILIZED MIXED FORMULATION: UNSTEADY BRINKMAN EQUATION 

The proposed formulation is based on the method of horizontal lines (also known as the Rothe 
method) |42j and variational multiscale formalism [p5]. To this end, we shall discretize the time 
interval of interest into iV instants denoted by t n , (n = 0, • • • , N). For simplicity, we shall assume 
uniform time steps, and shall denote the time step by At := t n — t n -\. However, it should be noted 
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that a straightforward modification can handle non-uniform time steps. We shall denote the time 
discretized version of a given quantity z(x,t) at the instant of time t n as follows: 



z {n) (x) « z(x,t n ) n = 0, 



N 



(6) 



In a semi-discrete formulation, the partial differential equation (which depends on both space 
and time) is first spatially discretized by (say) the finite element method. This results in a system 
of ordinary differential equations (ODEs) or a system of differential algebraic equations (DAEs), 
which are numerically integrated using appropriate time stepping schemes. On the other hand, in 
the method of horizontal lines (also known as the Rothe's method), the partial differential equation 
(which depends on both space and time) is temporally discretized using a time stepping scheme. 
This results in another system of partial differential equations, which depends only on spatial 
coordinates, and the finite element method or the finite difference method is typically employed 
to solve this resulting system of equations. Figure [T] gives a pictorial description of the method of 
horizontal lines. 



< 



< 



Solve equation ([7]) at time t = t n+ \ 



Solve equation (JTj) at time t = t n 



X 

Figure 1. Rothe's method: The governing equations are solved at each time level 
after discretizing the time interval into discrete time levels. 



Herein, we employ the backward Euler time stepping scheme, which is first-order accurate and 
unconditionally stable when applied to linear system of ordinary differential equations [48J. By 
employing the method of horizontal lines based on the backward Euler scheme to equation (jla|) - 

(lldp . the corresponding time discretized equations at time level t = t n+ \ can be written as follows: 
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«.(n+l) _ 7 ,(n) 

p — = -grad[p (n+1) ] - oV n+1 ) + div[/igrad[i; (ri+1) ]] + pb {n+l ^ in O (7a) 

div[^ n+1 )] =0 in (7b) 

„("+!) (a-) = w P(a; 5 t = t n+1 ) onT (7c) 

■u (0) (;e) = v (x) in ft (7d) 



Remark 3.1. Equation <\7ah can be rearranged to obtain the following equation: 

( At + ") V(n+1) = -S rad ^ (n+1) ] + div[/igrad[« (n+1) ]] + p (b^ + ^« (n) ) in Q (8) 

iime /eveZ t = t n+ i, the above equation ([8]) along with equations (|7bl )^( l7dj) resemble the well- 
known steady Brinkman equation with (modified) drag coefficient (a + -fa) and (modified) body 
force (V n+1 ) + ^u (n) ) . 



For convenience, let us define a, £(■, ■), and + \x) as follows: 

a := aAt + p (9a) 

£(w, q) := aw + At grad[q] — Aidivf/zgradftu]] (9b) 

b {n+l) :=Atb^+v^ (9c) 

Let L 2 (ft) denote the space of square integrable (scalar) functions defined on Q. The standard L 2 
inner product over domain K is denoted as (•, -)k- That is, 

(a;b) K = [ a-bdK Va,6eL 2 (ft) (10) 
Jk 

For simplicity, the subscript K will be dropped if the domain is whole of ft (that is, K = f2). In 
the remainder of the paper, we shall use the the following function spaces: 

V t := [v(x) | v(x) € (H\n)) nd , v(x) = v p (x, t) on r} (11a) 

W := [v(x) | G (H 1 (n)) nd , v(x) = on r} (lib) 

V := |p(a?) G L 2 (ft) | ^ p dft = o| (11c) 

Q:= |p(aj) G ^(O) | ^pdO = oj (lid) 

where iif 1 is a standard Sobolev space on ft [T7J. In the next subsection, we present a derivation 
of the proposed stabilization mixed formulation. 



Remark 3.2. If traction were prescribed on a part of the boundary, then pressure would be uniquely 
determined and we need not enforce J n p df2 = 0. However, in this paper, we do not consider the 
possibility of prescribing traction. 

3.1. A derivation of the proposed stabilized mixed formulation. We start with the classical 
mixed formulation for equations (|7a f )— (|7d|) . which can written as: Find v^ n+1 \x) £ Vt=t n+1 and 
pi n + l )(x) G V such that we have 

(w;av( n+ V) + At(grad[uj];/i grad[?/ n+1 )]) - At(div[u;];p( n+1 )) - At(q; div[v^ n+1 ^]) 
= {w;pb {n+1] ) Vio(a:) E W, q{x) £ V (12) 

As mentioned earlier, the classical mixed formulation has numerical instabilities, and hence the 
need for alternate (stabilized) formulations. To this end, let us divide the domain Q into "Nele" 
non-overlapping subdomains Q e (which, in the finite element context, will be elements) such that 

Nele 

n={Jti e (13) 

e=l 

The boundary of the element Vt e is denoted by T e := Cl e — £l e . We decompose the velocity field 
v( n+1 \x) into coarse-scale and fine-scale components, indicated as v c n+1 \x) and iA n (a;), re- 
spectively. That is, 

v (n+l) ( x j = v (n+l) ^ + v (n+l) ^ ^ 

Likewise, we decompose the weighting function w(x) into coarse-scale w c (x) and fine-scale Wf(x) 
components. At the outset, we recognize that this multi-scale decomposition can be done irre- 
spective of whether the problem is linear or non-linear. We further make an assumption that the 
fine-scale components vanish along each element boundary. That is, 

vf +1 \x) = w f (x) = on r e ; e = 1, • • • ,Nele (15) 

Let V c be the function space for the coarse-scale component of the velocity v c n+1 \x), and W c be 
the function space for w c (x). The spaces V c and W c are defined as 

V c := Vf, W c := W (16) 

where Vt and W are as defined earlier in equation (jllap and equation (jllbp . respectively. Let Vf be 
the function space for both the fine-scale component of the velocity v^ n+1 \x) and its corresponding 
weighting function Wf(x). It is defined as follows: 

V f := {v(x) E (H 1 (S})) nd | v(x) = on T e ; e = 1, • • • , Nele} (17) 
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Substituting equation (|M|) into equation (|T2|) . one obtains two subproblems. The coarse-scale and 
fine-scale subproblems are, respectively, defined as follows: Find v^ +l \x) G V c , v <y j l+1 \x) G V/, 
and p( n+1 \x) G V such that we have 

(w c ; av^ +1 ^ + (tu c ; d^. n+1) ) + At (grad[u> c ]; /i grad[^ +1 )]) 

+ At(grad[u; c ];// grad^ ( ™ +1) ]) - At (div[w c ];p( n+1 )) - At (g; div[i;( n+1 ); 



- At(q; div[vf +1) ]) = (w c ; pb (n+1) ) Vw c (x) G W c , q(x) G V (18a) 

(«;/; a«( n+1) ) + (to/; di>J. n+1) ) + At (grad[io/];/z grad[-y( n+1) ] 

+ At(grad[w / ];^grad[^ n+1) ]) - At(div[w/];p( n+1 )) = (to/; pb (n+1) ) Vto/(aj) G W/ (18b) 

The fine scale problem (|18bp can be solved (either approximately or exactly) and substituted into 
the coarse scale problem (|18a|) to eliminate the fine scale variables completely. But before we 
attempt to do that, it is necessary to interpolate the fine-scale variables. [We interpolate the 
coarse-scale variables using the standard Lagrangian shape functions.] The fine-scale variables are 
interpolated using bubble functions. A bubble function defined over an element is a function that 
vanishes on the boundary of the element. The fine-scale variables over an element are interpolated 
as follows: 

vf +l \x) = b e (x)(3, w f {x) = b e {x) 1 (19) 

where b e (x) is the bubble function while (3 and 7 are column vectors of size nd x 1. It will be 
easier to construct a bubble function on the reference element, and use the isoparametric mapping 
to define it with respect to the original coordinate system. Therefore, the bubble functions will be 
functions of £, which is the position vector in isoparametric coordinates. An illustration of all the 
different bubbles used in this paper is provided in Figure [2j 

Solving (|18bp for (3 and substituting into (|18ap . the proposed stabilized mixed formulation at 
time level t = i n +i we arrive at is the following (dropping the suffix c for convenience): Find 
v^ n+l \x) G Vt=t n+1 and p^ n+l \x) G Q such that we have 

(w;av( n+ V) + At(grad[™];^grad[<(/ n+1 )]) - At(div[w];p( n+1 )) - At(q; div[i/ n+1 )]) 

Nele Nele 

e=l e=l 

Vto(sB) G W, q(x) G Q (20) 
where the stabilization parameter t(x) is given by 

t(x) = b e {x) (^j b e (x) dnj (^j (fi At ||grad[6 e ]|j 2 + a{b e ) 2 ) dtoj (21) 
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where b e (x) is a bubble function. Thus we obtain a formulation which has additional stabilizing 
terms compared to classical mixed formulation. A systematic numerical implementation of the 
proposed formulation is outlined in Algorithm [TJ A few remarks about the stabilization parameter 
are in order. 

Remark 3.3. We have a > 0, fj, > 0, At > 0. Also, for x G Q e , {b e (x)) 2 > 0, ||grad[6 e ] || 2 > 0. 
The stabilization parameter t(x) is thus well-defined for all permissible a, p and At since the 
integrand (p At ||grad[6 e ] || 2 + a(b e ) 2 ) > 0. 

Remark 3.4. On defining 

1 



t{x) dQ (22) 



<> 



'avg •— /„ 

meas(ir 

one can get bounds on r avg for some special cases of interest. The Poincare inequality [39] implies 
that 

[ (b e (x)) 2 an < j(n e ) [ || g rad[6 e ]|| 2 dn (23) 

where 7(^ e ) is a domain- dependent positive constant. Note that, in the above equation, we have 
used the fact that the bubble function vanishes on the boundary. In particular, we have assumed 
that b e (x) € H^(Q e ). The elements in a typical finite element mesh will all be convex domains. For 
such domains, there are good estimates for the parameter 7 (for example, see Reference [50],). The 
Cauchy-Schwarz inequality on the other hand implies 

(^J b e (x)dt?J <meas(S! e )(y (b e (x)) 2 dttj (24) 

The inequalities (f23|) and (f24"|) immediately yield 

0<r avg < (a + n At 7 (0 e ))- 1 (25) 

Specifically, for the unsteady Darcy equation (fj, = 0) and the unsteady Stokes equation (a = p), we 
have 

< r avg < (a)- 1 if n = (26a) 
< r avg < (p + fi At 7(0 e ))" 1 if a = p (26b) 

Remark 3.5. Readers consulting the references provided herein are forewarned of a few errors and 

oversights in [51J . The boundary condition prescribed in equation (3) of [51J is appropriate to the 

Darcy equation and not to the Brinkman equation. A similar comment applies to equation (4) where 

the function space for the velocity is defined; such a space is valid for Darcy equation but not for 

the Brinkman equation. Moreover, in equation (4), it ought to read as trace(i;-n) = ip E H~ l / 2 {T). 
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Finally, the stabilization parameter defined in equation (30) will make sense only if we understand 
fiAv to mean pAv + /udiv[grad[t>] T ]. 



Algorithm 1 Implementation of the proposed formulation 
1: Input: Initial condition vq, Time period of integration T, time step At 
2: Set «W = v 
3: Set t = 
4: while t < T do 

5: At = min(At,T - t), t = t + At 

6: Using v^ n \ solve equation ([20]) to obtain v^ n+1 \ p( n+1 ) 

7: Set t>W = -y( n+1 ) 

8: end while 



Remark 3.6. in f/ie implementation of the algorithm, all the terms in equation (|20p need not be 
evaluated repeatedly at each time step since most of them do not depend on the temporal variable. 
In fact, only the terms involving b need to be evaluated repeatedly. 

4. SPECIAL CASE : A STABILIZED MIXED FORMULATION FOR UNSTEADY 

DARCY EQUATION 

We now present a stabilized mixed formulation for the unsteady Darcy equation following the 
same procedure as outlined in the previous section. Most of the quantities appearing have been 
defined previously. The ones appearing for the first time are defined here. The boundary T is 
divided into two parts, denoted by T v and T p , such that T v n T p = and T v U T p = T. F p is the 
part of the boundary on which the pressure p is prescribed, and T v is that part of the boundary 
on which the normal component of the velocity is prescribed. The prescribed normal component 
of the velocity on the boundary is denoted by ip(x,t), and the prescribed pressure is denoted by 
Po(x,t). The unsteady governing equations can be written as 

-grad[p] -av + pb in x (0, T) (27a) 

infix (0, T) (27b) 

--Po(x,t) onr p x(0,T) (27c) 

n(x) = i(>(x, t) on T v x (0, T) (27d) 

= vq(x) in Q, (27e) 

where n{x) is the unit outward normal on the boundary. 
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dv 

div[tt] = 
p(x,t) = 
v(x, t) ■ 
v(x,0) -- 



Remark 4.1. It is important to note that if T v = T, then for a given f2 (and hence given Y) ip 
cannot be arbitrarily specified. Using the divergence theorem it is apparent that ip must meet the 
following compatibility condition: 



tji dr = 



r=r „ 



(28) 



The relevant function spaces for the velocity and pressure fields in case of the unsteady Darcy 
equation are defined as follows: 



V t := {v(x) | v(x) G i? 1 (div;^),-j;(a;) ■ n(x) = ip(x,t) on T v } 
W := {v(x) | v(x) G iJ^div; n),v(x) ■ n(x) = on T v } 
V := {p(x) | p(x) G L 2 (Q)} , Q := {p(x) \ p(x) G H 1 ^)} 
where ff 1 (div;0) is defined by 

H^divjn) := |w(aj) G (L 2 (ft))" d | div[t>] G L 2 (Q)} 
Remark 4.2. If we have T v = T, then 

V := \p(x) G L 2 (9) | f pdfi = ol, Q := \p(x) G F X (0) | / p dO = 

For convenience, let us define the operator £.{■, •), a and b^ n+1 \x) to be 

£(io, g) := aw + At gradfg] 
a := p + aAt 



~(n+l) 



Atb («+1) + v (n) 



(29a) 
(29b) 
(29c) 

(30) 
(31) 

(32a) 
(32b) 
(32c) 



A stabilized mixed formulation for unsteady Darcy equation at time level t = t n +i can be written 
as: Find v( n+1 \x) G Vt=t n+1 and p^ n+l \x) G Q such that we have 

(w;av^ n+1 ^ - At(div[w];p(" +1 )) - At(q; dW[v ( - n+1 ^}) 

C{w,q)]T{x)a- 1 £( V ( n+1 \p( n+i y 



(33) 



(w;pb 



+ At(wn;pD (n+1) )n 

Vw(aj) G W, g(as) G Q 



£(«;, g); r(a;)o! 1 / o6^ l+1 ' ) 



where the stabilization parameter t{x) is given by 



t(x) = b e (x) 



b e (x) an 



ke\2 



(34) 



1 



For simplicity, one can use a constant representative value of t(x) = r avg = — for all elements. A 
careful discussion on the use of the representative value for T avg can be found in References 
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Remark 4.3. While this formulation requires the pressure p G H 1 ^), it is possible to use a classical 
mixed formulation with different interpolations for the velocity and pressure so that we only require 
p € L 2 (Q). But the relaxation in the requirement comes at the cost of increasing complexity and 
greater difficulty in implementation. Moreover, the advantage is short-lived because as we proceed 
in the hierarchy to a Brinkman model, we will require further smoothness from the pressure field in 
any case. 

5. REPRESENTATIVE NUMERICAL RESULTS 

In this section, we shall show the performance of the proposed stabilized formulation using 
representative numerical examples. We also perform spatial and temporal numerical convergence 
studies. In all the numerical simulations, we have employed the equal-order interpolation for the 
velocity and pressure. The drag coefficient, density, and viscosity are all taken asa = p = /i = l 
for all the cases. 



Part 1. UNSTEADY DARCY EQUATION 

5.1. One-dimensional problem. Consider the domain to be of unit length. The body force, 
boundary conditions, and initial conditions are, respectively, given by 

b(x, t) = - sin(t) + -Ax 3 sin(t) + cos(t) (35a) 
v(0,t) =v(l,t) =cos(i) (35b) 
v(x,0) = l (35c) 

For uniqueness, the pressure is prescribed at x = 0, and the value is taken to be zero. The analytical 
solution is given by 

v(x,t) = cos(i) (36a) 

p(x, t) = — x 4 sin(t) (36b) 

The above problem is solved using the proposed formulation, and the computational domain is di- 
vided into 20 equal-sized two- node linear finite elements. In Figures and H] the obtained numerical 
solution is compared with the exact solution for time instants at T = 1 and T = 2, respectively. 

5.1.1. Temporal numerical convergence studies. To study the temporal convergence, we compared 

the numerical solutions obtained at time T = 1 using successive time steps of magnitude 10 -2 , 
10~ 2 10~ 2 10~ 2 10" 2 

, , and with a mesh of 128 elements. The terminal rate of convergence is 

2 4 ' 8 16 

found to be approximately first order (see Figure [5|) , which is what is expected of a backward Euler 
time discretization scheme. 
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5.1.2. Spatial numerical convergence studies. To study the spatial convergence, we compared the 
numerical solutions obtained at time T = 1 using hierarchical meshes (consisting of 8, 16, 32, 64 
and 128 elements) with a time step of 10 -4 . The terminal rates of convergence are found to be 
(approximately) first order for the pressure and second order for the velocity, which are shown 
Figure 

5.2. Two-dimensional problem. The computational domain is a bi-unit square. The initial 
condition for the velocity is taken to be zero, and homogeneous boundary conditions for the normal 
component of the velocity are prescribed on the whole boundary. The test problem is pictorially 



defined in Figure [6l The body force is given by 

b x (x, y, t) = (2 exp(t) — 1) sin(27rx) sin(27ry) + sin(7rt)y(l — y)(l — 2x) (37a) 

b y (x,y,t) = (2exp(i) - 1) (cos(2vrx) cos(2vry) - cos(2vrx)) +sin(vrt)x(l - x)(l - 2y) (37b) 

The analytical solution is given by 

v x (x, y, t) = (exp(i) — 1) sin(27nr) sin(27ry), (38a) 

v y (x,y,t) = (exp(i) — 1) (cos(27rx) cos(27ry) — cos(27rx)) (38b) 

p{x, y, t) = x(l - x)y(l - y) sin(Trt) (38c) 



This problem was solved in turn on a T3, Q4 and an unstructured T3 grid respectively (see Figure 
[7J. A grid of size 31 x 31 was used for the case of T3 and Q4 elements while 448 elements were 
used in the case of the unstructured T3 grid. A comparison of the exact and numerical solution at 
time T = 0.5 is presented in the contour plots in Figures 181- 1 101 

5.2.1. Spatial convergence. The spatial convergence was studied by comparing numerical solutions 
obtained at time T = 0.5 using successively finer meshes with 5, 9, 17 and 33 nodes along each side 
respectively with a time step of magnitude 10~ 3 . The terminal rate of convergence is found to be 
approximately second order for both pressure and velocity with both four-node quadrilateral and 
three- node triangular finite elements (see Figure ITTj) . 

Part 2. UNSTEADY BRINKMAN EQUATION 

5.3. One-dimensional problem. Consider the domain to be of unit length. The body force, 
boundary conditions, and initial conditions are, respectively, given by 

b(x,t) = cos(107r£) + 67T exp(t) cos(67rx) — 107r sin(107rf) (39a) 
u(0, t) = v(l, t) = cos(lOvrt) (39b) 
u(s,0) = l (39c) 

15 



For uniqueness, the pressure is prescribed at x = 0, and the value is taken to be zero. The analytical 
solution is given by 



v(x,t) = cos(lOvrt) (40a) 
p(x, t) = exp(i) sin(67ra;) (40b) 

The above problem is solved using the proposed formulation, and the computational domain is 
divided into 40 equal-sized two-node linear finite elements. In Figures [T2] and [T3] the obtained 
numericalsolution is compared with the exact solution for time instants at T = 0.1 and T = 0.2, 
respectively. 



5.3.1. Temporal numerical convergence studies. To study the temporal convergence, we compared 

10~ 2 

the numerical solutions obtained at time T = 1 using successive time steps of magnitude — - — , 
10~ 2 10~ 2 10" 2 

, and with a mesh of 128 elements. The terminal rate of convergence is found to 

4 ' 8 16 B 

be approximately first order (see Figure [14"]) , which is what is expected of a backward Euler time 
discretization scheme. 

5.3.2. Spatial numerical convergence studies. To study the spatial convergence, we compared the 
numerical solutions obtained at time T = 0.01 using hierarchical meshes (consisting of 16, 32, 64, 
128 and 256 elements) with a time step of 10 -5 . The terminal rates of convergence are found to 
be (approximately) first order for the pressure and second order for the velocity, which are shown 
Figure [TU 

5.4. Two-dimensional problem. The computational domain is a bi-unit square. The initial 
condition for the velocity is taken to be zero, and velocity is prescribed everywhere on the boundary. 
The test problem is pictorially defined in Figure . The body force is given by 

b x (x, y, t) = (2 exp(t) — 1) sin(27ra;) sin(27ry) + 8ir 2 (exp(t) — 1) sin(27rx) sm(2iry) 

+ 7r exp(i) cos(7rx) sin(7ry) (41a) 

b y (x, y, t) = (2 exp(t) — 1) (cos(27rx) cos(27ry) — cos(27rx)) + 8tt 2 (exp(i) — 1) cos(27rx) cos(27ry) 

— 4tt 2 (exp(t) — 1) cos(27tj;) + 7rexp(t) cos(7ry) sin(7rx) (41b) 

The analytical solution is given by 

v x (x, y, t) = (exp(t) — 1) sin(27ra;) sin(27ry), (42a) 
v y (x,y,t) = (exp(i) — 1) (cos(27rx) cos(2-7ry) — cos(27rx)) (42b) 
p(x, y, t) = sin(7nr) sin(7ry) exp(i); (42c) 
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This problem was solved in turn on a (six-node triangle element) T6 and (nine-node quadrilateral 
element) Q9 grid respectively (see Figure^]). A grid of size 31 x 31 was used for the case of Q9 while 
a 41 x 41 grid was used in the case of the T6 element. A comparison of the exact and numerical 
solution at time T = 0.20 is presented in the contour plots in Figures [T31 and [TBI 

Remark 5.1. The T3 and Q4 elements perform satisfactorily for the unsteady Darcy equation. 
However in case of the unsteady Brinkman equation, they perform poorly with respect to computation 
of the pressure, mainly because the second order term (Laplacian) is identically zero. Hence we have 
used higher order elements, namely the T6 and Q9 elements for the unsteady Brinkman equation. 
However, with higher order elements, the choice of bubble functions is quite critical and is by itself 
an important question. In this case, with better choice of bubbles, the results should improve. 

5.4.1. Spatial convergence. The spatial convergence was studied at two different time instants, 
T = 0.01 and T = 0.1, by comparing numerical solutions obtained using successively finer meshes 
with a constant time step. For the case of T = 0.01, we use a time step At = 10 -4 while for 
T = 0.1, At = 10~ 3 . In case of the T6 element, meshes consisting of 5, 9, 17 and 33 nodes along 
each side were used in succession while for the Q9 element, in turn, it consisted of 3, 5, 9, and 17 
nodes along each side. The terminal rate of convergence is found to be approximately second order 
for both pressure and velocity with both Q9 and T6 finite elements (see Figure [T71 and fl8|) . 

6. CONCLUSIONS 

In this paper, we have presented a stabilized mixed formulation for unsteady Brinkman equation. 
The proposed formulation is consistently derived based on the variational multiscale formalism 
combined with the method of horizontal lines. The derivation does not require the assumption that 
the fine-scale terms do not depend on time, which is required under many of the existing stabilized 
variational multiscale formulations for transient problems. Under the proposed formulation, it has 
been shown numerically that the equal-order interpolation for the velocity and pressure is stable 
(which is not the case with the classical mixed formulation). The performance of the proposed 
formulation is illustrated using various representative numerical examples. Spatial and temporal 
numerical convergence studies are also performed, and the proposed formulation performed well. 

7. APPENDIX 

In Appendix, we present a finite element discretization of Av and Aw. For comparison and 
completeness, we also present the conventional way of deriving stabilized mixed formulation based 
on the variational multiscale approach. 
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7.1. Finite element discretization of Av and Aw terms. We first introduce some notation 
and definitions, which be useful in writing the "stiffness" matrices in a systematic manner. Let A 
and B be matrices of size n x m and p x q, respectively. 

" 61,1 • • • h, q 



A = 



0-1,1 ■ ■ ■ ai,m 



CL n \ . . . CL r , 



B = 



b p ,i ... b. 



The Kronecker product of these matrices is an np x mq matrix, and is defined as 

ai,iB . . . ai :Tn B 

AQ B := j •. : 

a n ,iB . . . a n , m B 

Note that the Kronecker product can defined for any two given matrices (irrespective of their 
dimensions). The vec[-] operator is defined as 

01,1 



vec[A] :- 



0-71,1 



Some relevant properties of Kronecker product and vec[-] operator are as follows. 

• vec[ACB] = (B T A) vec[C] 

• (A B) {C D) = {AC BD) 

• vec[A + B] = vec[A] + vec[S] 

For convenience and effective computer implementation, we shall represent a three-dimensional 
array as a two-dimensional matrix. Consider a three-dimensional array H of size m x n x p, and 
its corresponding matrix representation is denoted by mati(-) and is defined as follows: 

( H\u ■■■ Hi n i Hn p ■■■ Hi np ^ 

mati[if] := [ ' • . [ \ \ \ ' • . \ 

\H m ii ■ ■ ■ H mn i H m ip ■ ■ ■ H mn pJ 

Using the above notation, the operation Ui = Hij^Ajk (which is written in indicial notation) can 
be compactly written as 



u = mati[i3"]vec[A] 
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(43) 



Let C denote the position vector in the reference element. The row vector containing shape 
functions is denoted by N, which is a function of £. The derivatives and hessian of the vector of 
shape functions with respect to C are ) respectively, denoted as DN and DDN. That is, 

dNi d 2 N- 
(DN h= ^. and { DDN )ijk = ^ (44) 

The velocity vector is approximated as follows: 

v(x) = v T N T (C(x)) (45) 

where v denotes the matrix containing nodal velocities. (Note that nodal quantities are denoted 
with superposed hats.) Based on the isoparametric mapping, the Jacobian matrix and div[«7] can 
be written as 

J = x T (DN) (46) 
divfj- 1 ] = -J- 1 & T mati[DD.ZV]vec[J- 1 J- T ] (47) 

where x is a matrix containing nodal coordinates. For convenience, let us define 

y = (I — DNJ~ T x T )mati [DDN]vec[J~ 1 J~ T ] (48) 

Using the above notation, a finite element approximation for Av can be compactly written as 

Av = (y T Q I)vec[v T ] (49) 

Similarly, one can approximate Aw. 

7.2. Conventional way of deriving a stabilized formulation for transient Brinkman equa- 
tion. We define the relevant function spaces, following Hughes |48| . as follows (for a more rigorous 
definition of the function spaces, the reader may consult Reference |49j ) : 

V* := lv(-,t) | v(-,t) G (H 1 (n)) nd , v(x, t) = v p (x, t) on r} (50a) 



V* := jp(-,t) G L 2 (tt) j y"pdfi = o| 
Q t := |p(-,t) G ff x (0) | J pdn = 0^ 



(50b) 
(50c) 



We again start with the classical mixed formulation for the governing equations: Find v(x,t) G V* 
and p(x, t) G V 1 such that we have 

dv \ 

w 'iP~H7 + (w;av) + (grad[u>]; /igrad[u]) - (div[io];p) - (q;div[v}) 



dt 



{w;pb) Vto(aj) G W, q(x) G T 3 (51) 
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We decompose the velocity as follows: 



v(x, t) = v c (x, t) + Vf(x) 



(52) 



As pointed out earlier in Section [TJ the above assumption that Vf(x) is independent of time, while 
mathematically plausible, is philosophically undesirable and does not seem to have a physical basis. 
Substituting (|52p into (|5ip . we obtain the coarse-scale and fine-scale subproblems, which can be 
written as follows: Find v c (x,t) G V*, Vf(x) G V/ and p(x,t) G V 1 such that 



w c ;p- 



dv c 
dt 



(w c ;av c ) + (w c ;av f ) + (grad[w c ]; p grad[v c ]) + (grad[w c ]; p grad[v/]) 



(div[tu c ];p) - (q;div[v c ]) - (g; div[t>/]) = (w c ;pb) Mw c (x) G W c , g(sc) G T 3 (53a) 



w f'P-Qf j + ('"'/; at, c) + (wf,avf) + (grad[«7/];// grad[« c ]) + (grad[w/]; grad[v/]) 
(div[w f ];p) = (w f ;pb) Vwf(x) E W/ (53b) 



The remaining part of the derivation is as follows: we interpolate the fine-scale variables (vf(x) 
and Wf(x)) in terms of the bubble function b e (x). We then solve the fine-scale problem to obtain 
the fine-scale component of the velocity in terms of the coarse-scale component of the velocity and 
pressure. The obtained expression for the fine-scale component of the velocity is then substituted 
in the coarse-scale subproblem (|53ap to eliminate the fine-scale degrees-of-freedom and obtain a 
stabilized mixed formulation. For convenience, dropping the suffix c, the final form of the stabilized 
mixed formulation reads: Find v(x,t) G V* and p(x,t) G Q l such that we have 

dv \ 

w;p— + (w;av) + (grad[u/|; p grad[vj) - (div[w];p) - (q;div[v]) 



Nele 



dv 



( aw + gradfg] — //div[grad[«j]; f(x)p- 

Nele 

(aw + grad[q] — /xdiv[grad[t«]; t(x) (av + grad[p] — /udiv[grad[i;])) $ 

e=l 

Nele 

(w; pb) — (aw + grad[g] — /xdiv[grad[t«]; f(x)pb) ne Vw(x) G W, q(x) G Q 



Hi'- 



(54) 



e=l 



where the stabilization parameter in this case is defined as follows: 



f(x) = b e (x) 



b e (x) dn 



20 



(p ||grad[6 e ]|| 2 + a(6 e ) 2 ) 



n -l 



(55) 



A semi-discretization of equation (|54j) yields a system of ordinary differential equations (ODEs) 
of the following general form [48]: 



M vv 


M vp " 




V 


+ 


K vv 


K vp " 




V 




V 


















M pv 







_P_ 




K pv 


Kp P . 




_P_ 




fp. 



(56) 



The above system of ODEs can be solved by employing any of the known time stepping schemes 
(e.g., Runge-Kutta methods, generalized trapezoidal family of time stepping schemes). 
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Figure 2. This figure shows the bubble functions for various finite elements. The 
coordinate axes of the reference/master element are denoted by £i and £2- The 
locations of the degrees-of-freedom for velocity and pressure nodes in each element 
are also indicated. 
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Figure 3. One-dimensional test problem (unsteady Darcy equation): In this figure 
the obtained numerical solution is compared with the analytical solution for the 
time instant at T = 1. The computational mesh consisted of 20 equal-sized two- 
node linear finite elements. 
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Figure 4. One-dimensional test problem (unsteady Darcy equation): In this figure 
the obtained numerical solution is compared with the analytical solution for the 
time instant at T = 2. The computational mesh consisted of 20 equal-sized two- 
node linear finite elements. 
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(a) Temporal convergence using 128 elements after time T = 1 
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(b) Spatial convergence using At = 10 4 after time T — 1 

Figure 5. One-dimensional test problem (unsteady Darcy equation): Temporal 
and spatial convergence. Note that the rate of convergence is in good agreement 
with the theory. 
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Figure 6. Two-dimensional problem : A pictorial description of the test problem. 
The computational domain is a bi-unit square. For unsteady Darcy equation: Ho- 
mogeneous normal components of the velocity are prescribed on the whole boundary 
(that is, 4>{x,t) = 0). The components of the velocity vector are, respectively, de- 
noted by v x and v y so that v y = on Ti, T3 while v x = on T2, For unsteady 
Brinkman equation: Velocity is prescribed completely on the boundary, so that, 
v x = v y = on Ti, T3 while v x = 0, v y = (exp(£) — 1) (cos(27ry) — 1) on T2, T4. To 
ensure uniqueness of the solution, we enforce pressure to be zero at the origin O in 
both cases. 
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Figure 7. Two-dimensional test problem: Computational meshes used in the nu- 
merical simulations. 
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Figure 8. Two-dimensional test problem (unsteady Darcy equation) with three- 
node triangular finite elements: Comparison of analytical and numerical solutions 
at final time T = 0.5. An equally spaced 31 x 31 grid was used with At = 10 -3 . 
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Figure 9. Two-dimensional test problem (unsteady Darcy equation) with four- 
node quadrilateral finite elements: Comparison of analytical and numerical solutions 
at final time T = 0.5. An equally spaced 31 x 31 grid was used with At = 10 -3 . 
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Figure 10. Two-dimensional test problem (unsteady Darcy equation) with un- 
structured grid of three-node finite elements: Comparison of analytical and numeri- 
cal solutions at final time T = 0.5. 1792 elements were used with At = 10 -3 . 
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Figure 11. Two-dimensional test problem (unsteady Darcy equation): The figure 
presents spatial convergence of the proposed formulation using four-node quadrilat- 
eral elements (top) and three- node triangular elements (bottom). The time step is 
taken to be At = 10 -3 , and the total time is taken to be T = 0.5. 
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FIGURE 12. One-dimensional test problem (unsteady Brinkman equation): In this 
figure the obtained numerical solution is compared with the analytical solution for 
the time instant at T = 0.1. The computational mesh consisted of 40 equal-sized 
two-node linear finite elements. 
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FIGURE 13. One-dimensional test problem (unsteady Brinkman equation): In this 
figure the obtained numerical solution is compared with the analytical solution for 
the time instant at T = 0.2. The computational mesh consisted of 40 equal-sized 
two-node linear finite elements. 
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(a) Temporal convergence using 128 elements after time T = 1 
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(b) Spatial convergence using At =10 after time T = 0.01 

Figure 14. One-dimensional test problem (unsteady Brinkman equation): Tempo- 
ral and spatial convergence. Note that the rate of convergence is in good agreement 
with the theory. 
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Figure 15. Two-dimensional test problem (unsteady Brinkman equation) with six- 
node triangular finite elements: Comparison of analytical and numerical solutions 
at final time T = 0.2. An equally spaced 41 x 41 grid was used with At = 10 -3 . 
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Figure 16. Two-dimensional test problem (unsteady Brinkman equation) with 
nine-node quadrilateral finite elements: Comparison of analytical and numerical 
solutions at final time T = 0.2. An equally spaced 31 x 31 grid was used with 
At = 1(T 3 . 
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Figure 17. Two-dimensional test problem (unsteady Brinkman equation): The 
figure presents spatial convergence of the proposed formulation using six-node tri- 
angular finite elements. The time step and the total time are taken to be At = 10~ 4 
and T = 0.01 (top) and At = 10~ 3 and T = 0.10 (bottom), respectively. 
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Figure 18. Two-dimensional test problem (unsteady Brinkman equation): The fig- 
ure presents spatial convergence of the proposed formulation using nine-node quadri- 
lateral finite elements. The time step and the total time are taken to be At = 10~ 4 
and T = 0.01 (top) and At = 10 -3 and T = 0.10 (bottom), respectively . 
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